Forecasting Swiss Energy Demand

---

The Project

---

TBD

---

Data Procurement and Cleaning

---

Firstly, I start by loading the data. Looking at the first 5 rows:

glimpse(data)

## Rows: 461,471
## Columns: 2
## $ date <dttm> 2009-01-01 00:15:00, 2009-01-01 00:30:00, 2009-01-01 00:45:00, 2…
## $ kwh  <dbl> 1829966, 1715173, 1732582, 1706689, 1722373, 1691958, 1686226, 16…

There are no missing values. However, the frequency is too granular. For modelling purposes, I want to aggregate to daily values and filter out everything after 2019 to exclude the effect of the pandemic.

data <- data %>% 
  mutate(year = year(date),
         day_in_year = yday(date)) %>% 
  filter(year <= 2021) %>% 
  group_by(year, day_in_year) %>% 
  summarise(date = last(date) %>% as.Date(),
            mwh = sum(kwh)/1e3) %>% 
  ungroup()

summary(data)

##       year       day_in_year         date                 mwh        
##  Min.   :2009   Min.   :  1.0   Min.   :2009-01-01   Min.   :102894  
##  1st Qu.:2012   1st Qu.: 92.0   1st Qu.:2012-04-01   1st Qu.:141458  
##  Median :2015   Median :183.0   Median :2015-07-02   Median :153739  
##  Mean   :2015   Mean   :183.1   Mean   :2015-07-02   Mean   :155591  
##  3rd Qu.:2018   3rd Qu.:274.0   3rd Qu.:2018-10-01   3rd Qu.:171736  
##  Max.   :2021   Max.   :366.0   Max.   :2021-12-31   Max.   :222496

Now we have got ten years worth of training data (2009-2018) and one year worth of testing data (2019).

Temperature Data

Got this document from opendata.swiss.

This link led to this document, detailing various weather stations:

read_csv("Data/liste-download-nbcn-d.csv")

Using the URLs to access the actual weather data for selected weather stations:

temp_data <- read_csv("Data/liste-download-nbcn-d.csv") %>%
  filter(`station/location` %in% c("SMA", "GVE", "BAS", "DAV","LUG", "STG",
                                   "ALT", "NEU", "BER", "LUZ")) %>%
  mutate(station_data = map(`URL Previous years (verified data)`,
                            ~ read.csv(.x, sep = ";")),
         station_data = map(station_data,
                            ~ select(.x, date, tre200d0)),
         station_data = map(station_data,
                            ~ na_if(.x, "-")),
         station_data = map(station_data,
                            ~ .x %>%
                              mutate(date = as.character(date),
                                     tre200d0 = as.numeric(tre200d0)))) %>%
  select(Station, `station/location`, Canton, station_data) %>%
  unnest(station_data) %>%
  mutate(date = ymd(date)) %>%
  select(-Station, -Canton) %>%
  pivot_wider(names_from = `station/location`, values_from = "tre200d0")

Left-join temperature data to the clean electricity usage data and create average temperature column:

data <- data %>% 
  left_join(temp_data, by = "date") %>% 
  mutate(avg_temp = rowMeans(select(., ALT:STG)))

summary(data)

##       year       day_in_year         date                 mwh        
##  Min.   :2009   Min.   :  1.0   Min.   :2009-01-01   Min.   :102894  
##  1st Qu.:2012   1st Qu.: 92.0   1st Qu.:2012-04-01   1st Qu.:141458  
##  Median :2015   Median :183.0   Median :2015-07-02   Median :153739  
##  Mean   :2015   Mean   :183.1   Mean   :2015-07-02   Mean   :155591  
##  3rd Qu.:2018   3rd Qu.:274.0   3rd Qu.:2018-10-01   3rd Qu.:171736  
##  Max.   :2021   Max.   :366.0   Max.   :2021-12-31   Max.   :222496  
##       ALT              BAS              BER               DAV         
##  Min.   :-10.30   Min.   :-12.40   Min.   :-14.200   Min.   :-22.200  
##  1st Qu.:  4.50   1st Qu.:  5.50   1st Qu.:  3.500   1st Qu.: -1.300  
##  Median : 10.80   Median : 11.40   Median :  9.700   Median :  4.500  
##  Mean   : 10.54   Mean   : 11.28   Mean   :  9.661   Mean   :  4.268  
##  3rd Qu.: 16.50   3rd Qu.: 17.10   3rd Qu.: 15.800   3rd Qu.: 10.200  
##  Max.   : 28.50   Max.   : 29.20   Max.   : 27.700   Max.   : 21.400  
##       GVE            LUG             LUZ              NEU        
##  Min.   :-9.2   Min.   :-5.30   Min.   :-11.60   Min.   :-10.20  
##  1st Qu.: 5.3   1st Qu.: 7.10   1st Qu.:  4.30   1st Qu.:  5.00  
##  Median :11.1   Median :13.20   Median : 10.60   Median : 10.90  
##  Mean   :11.3   Mean   :13.37   Mean   : 10.44   Mean   : 10.99  
##  3rd Qu.:17.3   3rd Qu.:19.60   3rd Qu.: 16.40   3rd Qu.: 17.00  
##  Max.   :29.5   Max.   :28.80   Max.   : 28.00   Max.   : 28.80  
##       SMA              STG             avg_temp      
##  Min.   :-13.50   Min.   :-14.900   Min.   :-12.060  
##  1st Qu.:  4.10   1st Qu.:  3.300   1st Qu.:  4.098  
##  Median : 10.30   Median :  9.400   Median : 10.155  
##  Mean   : 10.13   Mean   :  9.051   Mean   : 10.104  
##  3rd Qu.: 16.10   3rd Qu.: 14.900   3rd Qu.: 16.040  
##  Max.   : 27.80   Max.   : 27.100   Max.   : 26.970

Data is clean and ready to be used for model training. I split it into pre- and post-COVID for later test of the model regarding robustness during corona years:

data <- data %>% 
  filter(year <= 2019)

data_corona <- data %>% 
  filter(year > 2019)

Only the data prior to 2019 is used for the modelling.

---

Exploratory Data Analysis

---

data %>% 
  mutate(day_in_week = wday(date),
         day = wday(date, label = T),
         weekend = case_when(
           day_in_week %in% c(7) ~ "Saturday",
           day_in_week %in% c(1) ~ "Sunday",
           TRUE ~ "Work Day") %>% 
           factor(levels = c("Work Day", "Saturday", "Sunday"))) %>% 
  ggplot(aes(day_in_year, mwh, colour = weekend)) +
  geom_point(alpha = 0.4) +
  labs(title = "On Sundays, You Shall Not Work: Swiss End-User Electricity Consumption",
       subtitle = "For each day in the year (1-365), 11 daily consumption points from the period 2009 to 2019 are plotted.\nColouring indicates if a day lies on a weekend or not.",
       caption = "Data Source: Swissgrid, Period 2009-2019",
       x = "Day Of The Year (1-365)",
       y = "Daily Consumption",
       colour = NULL) +
  scale_x_continuous(breaks = c(1, seq(50,366,50))) +
  scale_y_continuous(labels = comma_format(suffix = " MWh")) +
  ggsci::scale_colour_jama() +
  theme_light() +
  theme(plot.title = element_text(face = "bold", size = 14),
        plot.subtitle = element_text(face = "italic", size = 10,
                                     colour = "grey50"),
        plot.caption = element_text(face = "italic", size = 9,
                                    colour = "grey50"),
        legend.position = "right")

From the above chart, it becomes clear that the demand follows a strong intra-year pattern. Additionally, there is a seasonality component contained within each week. Monthly patterns are not discernible.

data %>% 
  filter(year == 2010) %>%
  ggplot(aes(date, mwh)) +
  geom_line() +
  labs(title = "Daily Energy Consumption in 2010",
       y = "Daily Consumption",
       x = NULL) +
  scale_y_continuous(labels = comma_format(suffix = " MWh")) +
  theme_light() +
  theme(plot.title = element_text(face = "bold", size = 14),
        plot.subtitle = element_text(face = "italic", size = 10,
                                     colour = "grey50"),
        plot.caption = element_text(face = "italic", size = 9,
                                    colour = "grey50"),
        legend.position = "right")

Temperature development at different stations:

data %>% 
  filter(year == 2019) %>% 
  select(-c(year, day_in_year, mwh)) %>% 
  pivot_longer(-date) %>% 
  ggplot(aes(date, value, colour = name)) +
  geom_line(size = 0.75, alpha = 0.5) +
  theme_bw()

Average temperature development and electricity usage:

data %>% 
  select(date, avg_temp, mwh) %>%  
  mutate(day_in_week = wday(date),
         day = wday(date, label = T),
         weekend = case_when(
           day_in_week %in% c(7) ~ "Saturday",
           day_in_week %in% c(1) ~ "Sunday",
           TRUE ~ "Work Day") %>% 
           factor(levels = c("Work Day", "Saturday", "Sunday"))) %>% 
  ggplot(aes(avg_temp, mwh, colour = weekend)) +
  geom_point(alpha = 0.5, size = 1) +
  geom_smooth(size = 0.5, se = F) +
  labs(title = "Relation Between Energy Usage And Temperature",
       subtitle = "TBD",
       y = NULL,
       x = "Average Temperature",
       colour = NULL) +
  scale_y_continuous(labels = scales::comma_format(suffix = " MWh")) +
  scale_x_continuous(labels = scales::comma_format(suffix = "°")) +
  ggsci::scale_colour_jama() +
  theme_bw() +
  theme(plot.title = element_text(face = "bold", size = 14),
        plot.subtitle = element_text(face = "italic", size = 10,
                                     colour = "grey50"),
        plot.caption = element_text(face = "italic", size = 9,
                                    colour = "grey50"),
        legend.position = "right")

---

Building A Model

---

In this post, I am not going to use usual time series forecasting techniques, but I’ll rely on machine learning methods instead.

Creating splits:

dt_split <- data %>% 
  time_series_split(date_var = date, assess = "1 year", cumulative = T)

dt_train <- training(dt_split)
dt_test <- testing(dt_split)

folds <- vfold_cv(dt_train, v = 5)

Visualising the train/test split, it becomes very clear, what the goal will be: Forecasting the year 2019, for which we have the true data and are able to compute metrics to evaluate model performance.

bind_rows(
  dt_train %>% mutate(id = "training"),
  dt_test %>% mutate(id = "testing")
) %>% 
  ggplot(aes(date, mwh, colour = id)) +
  geom_line() +
  scale_colour_manual(values = c("firebrick", "dodgerblue")) +
  labs(title = "Training/Testing Split") +
  scale_y_continuous(labels = comma_format()) +
  theme_light() +
  theme(plot.title = element_text(face = "bold", size = 14),
        plot.subtitle = element_text(face = "italic", size = 10,
                                     colour = "grey50"),
        plot.caption = element_text(face = "italic", size = 9,
                                    colour = "grey50"),
        legend.position = "right")

Creating a recipe for data preprocessing and squeezing as much information from the date column as possible.

model_rec <- recipe(mwh ~ .,
                    data = dt_train) %>%
  step_date(date) %>%
  step_holiday(date, holidays = timeDate::listHolidays("CH")) %>%
  step_mutate(month = lubridate::month(date),
              year_half = lubridate::semester(date) %>% as.factor,
              week_day = lubridate::wday(date),
              week_in_year = lubridate::week(date),
              quarter = lubridate::quarter(date),
              date_christmas = ifelse(between(day_in_year, 358, 366),
                                      1, 0),
              day_in_year = as.factor(day_in_year)) %>% 
  step_rm(date) %>% 
  step_dummy(all_nominal_predictors(), one_hot = TRUE) %>% 
  step_zv()

model_rec %>% prep() %>% juice() %>% glimpse()

## Rows: 3,652
## Columns: 411
## $ year                    <dbl> 2009, 2009, 2009, 2009, 2009, 2009, 2009, 2009…
## $ ALT                     <dbl> 0.7, 0.1, -1.3, -3.3, -3.0, -2.0, -1.6, -2.2, …
## $ BAS                     <dbl> 1.5, -1.0, -3.1, -3.7, -4.1, -4.2, -5.7, -6.4,…
## $ BER                     <dbl> -1.4, -2.3, -4.9, -6.0, -5.9, -5.0, -5.2, -6.1…
## $ DAV                     <dbl> -6.2, -9.7, -12.5, -8.9, -11.0, -10.3, -9.2, -…
## $ GVE                     <dbl> 0.5, 0.0, -2.2, -3.5, -4.3, -2.8, -2.9, -3.5, …
## $ LUG                     <dbl> 1.4, 1.4, 1.0, -0.2, 0.7, 0.1, 1.3, 2.5, 1.9, …
## $ LUZ                     <dbl> 0.5, -0.8, -3.6, -4.2, -5.1, -3.5, -3.9, -4.6,…
## $ NEU                     <dbl> -0.5, -1.3, -3.9, -5.1, -4.2, -3.6, -4.2, -5.1…
## $ SMA                     <dbl> -0.5, -2.2, -5.1, -5.3, -5.8, -4.6, -5.4, -6.5…
## $ STG                     <dbl> -1.2, -3.8, -6.1, -6.6, -8.1, -5.8, -6.7, -7.8…
## $ avg_temp                <dbl> -0.52, -1.96, -4.17, -4.68, -5.08, -4.17, -4.3…
## $ mwh                     <dbl> 143210.0, 154168.8, 160934.3, 154669.9, 187821…
## $ date_year               <dbl> 2009, 2009, 2009, 2009, 2009, 2009, 2009, 2009…
## $ date_CHAscension        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ date_CHBerchtoldsDay    <dbl> 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ date_CHConfederationDay <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ date_CHKnabenschiessen  <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ date_CHSechselaeuten    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ month                   <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ week_day                <dbl> 5, 6, 7, 1, 2, 3, 4, 5, 6, 7, 1, 2, 3, 4, 5, 6…
## $ week_in_year            <dbl> 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 3, 3…
## $ quarter                 <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ date_christmas          <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X1          <dbl> 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X2          <dbl> 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X3          <dbl> 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X4          <dbl> 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X5          <dbl> 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X6          <dbl> 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X7          <dbl> 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X8          <dbl> 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X9          <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X10         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X11         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0…
## $ day_in_year_X12         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0…
## $ day_in_year_X13         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0…
## $ day_in_year_X14         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0…
## $ day_in_year_X15         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0…
## $ day_in_year_X16         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1…
## $ day_in_year_X17         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X18         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X19         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X20         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X21         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X22         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X23         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X24         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X25         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X26         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X27         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X28         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X29         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X30         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X31         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X32         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X33         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X34         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X35         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X36         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X37         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X38         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X39         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X40         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X41         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X42         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X43         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X44         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X45         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X46         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X47         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X48         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X49         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X50         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X51         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X52         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X53         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X54         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X55         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X56         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X57         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X58         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X59         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X60         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X61         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X62         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X63         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X64         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X65         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X66         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X67         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X68         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X69         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X70         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X71         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X72         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X73         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X74         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X75         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X76         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X77         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X78         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X79         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X80         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X81         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X82         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X83         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X84         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X85         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X86         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X87         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X88         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X89         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X90         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X91         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X92         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X93         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X94         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X95         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X96         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X97         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X98         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X99         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X100        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X101        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X102        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X103        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X104        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X105        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X106        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X107        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X108        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X109        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X110        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X111        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X112        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X113        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X114        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X115        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X116        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X117        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X118        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X119        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X120        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X121        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X122        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X123        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X124        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X125        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X126        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X127        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X128        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X129        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X130        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X131        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X132        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X133        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X134        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X135        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X136        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X137        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X138        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X139        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X140        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X141        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X142        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X143        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X144        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X145        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X146        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X147        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X148        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X149        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X150        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X151        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X152        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X153        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X154        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X155        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X156        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X157        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X158        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X159        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X160        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X161        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X162        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X163        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X164        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X165        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X166        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X167        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X168        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X169        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X170        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X171        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X172        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X173        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X174        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X175        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X176        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X177        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X178        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X179        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X180        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X181        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X182        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X183        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X184        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X185        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X186        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X187        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X188        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X189        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X190        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X191        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X192        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X193        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X194        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X195        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X196        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X197        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X198        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X199        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X200        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X201        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X202        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X203        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X204        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X205        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X206        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X207        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X208        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X209        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X210        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X211        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X212        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X213        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X214        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X215        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X216        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X217        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X218        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X219        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X220        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X221        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X222        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X223        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X224        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X225        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X226        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X227        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X228        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X229        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X230        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X231        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X232        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X233        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X234        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X235        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X236        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X237        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X238        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X239        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X240        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X241        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X242        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X243        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X244        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X245        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X246        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X247        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X248        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X249        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X250        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X251        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X252        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X253        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X254        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X255        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X256        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X257        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X258        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X259        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X260        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X261        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X262        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X263        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X264        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X265        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X266        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X267        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X268        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X269        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X270        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X271        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X272        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X273        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X274        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X275        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X276        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X277        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X278        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X279        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X280        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X281        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X282        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X283        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X284        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X285        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X286        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X287        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X288        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X289        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X290        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X291        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X292        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X293        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X294        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X295        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X296        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X297        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X298        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X299        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X300        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X301        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X302        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X303        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X304        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X305        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X306        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X307        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X308        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X309        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X310        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X311        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X312        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X313        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X314        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X315        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X316        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X317        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X318        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X319        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X320        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X321        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X322        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X323        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X324        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X325        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X326        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X327        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X328        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X329        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X330        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X331        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X332        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X333        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X334        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X335        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X336        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X337        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X338        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X339        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X340        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X341        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X342        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X343        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X344        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X345        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X346        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X347        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X348        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X349        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X350        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X351        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X352        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X353        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X354        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X355        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X356        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X357        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X358        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X359        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X360        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X361        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X362        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X363        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X364        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X365        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ day_in_year_X366        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ date_dow_Sun            <dbl> 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0…
## $ date_dow_Mon            <dbl> 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0…
## $ date_dow_Tue            <dbl> 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0…
## $ date_dow_Wed            <dbl> 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0…
## $ date_dow_Thu            <dbl> 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0…
## $ date_dow_Fri            <dbl> 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1…
## $ date_dow_Sat            <dbl> 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0…
## $ date_month_Jan          <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ date_month_Feb          <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ date_month_Mar          <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ date_month_Apr          <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ date_month_May          <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ date_month_Jun          <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ date_month_Jul          <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ date_month_Aug          <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ date_month_Sep          <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ date_month_Oct          <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ date_month_Nov          <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ date_month_Dec          <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ year_half_X1            <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ year_half_X2            <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…

en_spec <- linear_reg(penalty = 0.02, mixture = 0.5) %>% 
  set_engine("glmnet") %>% 
  set_mode("regression")

xg_spec <- boost_tree(trees = 1000) %>% 
  set_engine("xgboost") %>% 
  set_mode("regression")

rf_spec <- rand_forest(trees = 500) %>% 
  set_engine("ranger") %>% 
  set_mode("regression")

Without using resamples or hyperparameter tuning, the workflows can immediately be entirely fit to the training data.

en_wf_fit <- workflow() %>% 
  add_recipe(model_rec) %>% 
  add_model(en_spec) %>% 
  fit(dt_train)

xg_wf_fit <- workflow() %>% 
  add_recipe(model_rec) %>% 
  add_model(xg_spec) %>% 
  fit(dt_train)

rf_wf_fit <- workflow() %>% 
  add_recipe(model_rec) %>% 
  add_model(rf_spec) %>% 
  fit(dt_train)

---

Evaluating Model Performance

---

With the fitted models, I can now make predictions on the holdout data:

models <- tibble(
  type = c("RF", "EN", "GB"),
  wflow = list(rf_wf_fit, en_wf_fit, xg_wf_fit),
  predictions = map(.x = wflow, .f = ~ augment(.x, dt_test))
)

models

From this, let’s extract the predictions for each model:

predictions <- models %>% 
  select(-wflow) %>% 
  unnest(predictions)

predictions

Defining a set of metric, I can now compare the performance of the three models with default parameters:

evaluation_metrics <- metric_set(rsq, rmse, mae)

predictions %>% 
  group_by(type) %>% 
  evaluation_metrics(truth = mwh, estimate = .pred) %>% 
  ggplot(aes(type, .estimate, colour = type)) +
  geom_point(size = 3, shape = 18) +
  facet_wrap(~ .metric, scales = "free_y") +
  theme_light() +
  theme(legend.position = "none",
        panel.grid.major.x = element_blank())

predictions %>% 
  group_by(type) %>% 
  evaluation_metrics(truth = mwh, estimate = .pred) %>% 
  select(-.estimator) %>% 
  pivot_wider(names_from = ".metric", values_from = ".estimate") %>% 
  arrange(-rsq)

Clearly, random forest and gradient boosting have outperformed elastic net.

predictions %>% 
  ggplot(aes(mwh, .pred)) +
  geom_point(alpha = 0.2, colour = "midnightblue", size = 2) +
  facet_wrap(~ type) +
  geom_abline(lty = "dashed", colour = "grey50") +
  scale_x_continuous(labels = scales::comma_format()) +
  scale_y_continuous(labels = scales::comma_format()) +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))

predictions %>% 
  rename(actual = mwh, prediction = .pred) %>% 
  select(type, date, actual, prediction) %>% 
  pivot_longer(-c(type, date)) %>% 
  ggplot(aes(date, value, colour = name)) +
  geom_line(size = 0.5) +
  labs(x = NULL,
       y = "mwh",
       colour = NULL) +
  facet_wrap(~ type, ncol = 1) +
  ggsci::scale_colour_jama() +
  scale_y_continuous(labels = scales::comma_format()) +
  theme_bw()

---

Doubling Down On Tree Ensembles

---

With the knowledge from above, let’s now tune hyperparameters for a random forest and gradient boosting model, to see whether we can fill the gaps of the standard models.

rf_spec <- rand_forest(mtry = tune(),
                       min_n = tune(),
                       trees = tune()) %>%
  set_engine("ranger") %>%
  set_mode("regression")

rf_wflow <- workflow() %>% 
  add_recipe(model_rec) %>% 
  add_model(rf_spec)

rf_grid <- grid_latin_hypercube(finalize(mtry(), dt_train),
                                min_n(),
                                trees(),
                                size = 50)

gb_spec <- boost_tree(mtry = tune(),
                      trees = tune(),
                      min_n = tune(),
                      tree_depth = tune(),
                      learn_rate = tune(),
                      loss_reduction = tune()) %>%
  set_engine("xgboost") %>%
  set_mode("regression")

gb_wflow <- workflow() %>% 
  add_recipe(model_rec) %>% 
  add_model(gb_spec)

gb_grid <- grid_latin_hypercube(finalize(mtry(), dt_train),
                                trees(),
                                min_n(),
                                tree_depth(),
                                learn_rate(),
                                loss_reduction(),
                                size = 100)

# Random Forest
unregister_dopar <- function() {
  env <- foreach:::.foreachGlobals
  rm(list=ls(name=env), pos=env)
}

start_time <- Sys.time()

cl <- makePSOCKcluster(6)
registerDoParallel(cl)

rf_tune <- tune_grid(object = rf_wflow,
                     grid = rf_grid,
                     resamples = folds)

stopCluster(cl)
unregister_dopar()

end_time <- Sys.time()
end_time - start_time

# Gradient Boosting
unregister_dopar <- function() {
  env <- foreach:::.foreachGlobals
  rm(list=ls(name=env), pos=env)
  }

start_time <- Sys.time()

cl <- makePSOCKcluster(6)
registerDoParallel(cl)

gb_tune <- tune_grid(object = gb_wflow,
                     grid = gb_grid,
                     resamples = folds)

stopCluster(cl)
unregister_dopar()

end_time <- Sys.time()
end_time - start_time

Comparing the model performances of both random forest and gradient boosting:

rf_tune %>% 
  show_best(metric = "rsq")

gb_tune %>% 
  show_best(metric = "rsq")

Gradient boosting performed slightly better. Therefore, it will be the final model.

Fitting onto the entire training data set:

gb_final_fit <- gb_wflow %>% 
  finalize_workflow(select_best(gb_tune, metric = "rsq")) %>% 
  last_fit(dt_split)

Evaluating out-of-sample performance:

gb_final_fit %>% 
  collect_predictions() %>% 
  ggplot(aes(mwh, .pred)) +
  geom_point(alpha = 0.5, size = 1.5, colour = "midnightblue") +
  labs(title = "Out-Of-Sample Gradient Boosting Model Performance",
       y = "Estimate",
       x = "Truth") +
  geom_abline(colour = "grey50", lty = "dashed") +
  scale_y_continuous(labels = scales::comma_format(suffix = " MWh")) +
  scale_x_continuous(labels = scales::comma_format(suffix = " MWh")) +
  theme_bw() +
  theme(plot.title = element_text(size = 14, face = "bold"),
        plot.subtitle = element_text(face = "italic", size = 12,
                                     colour = "grey50"))

Out-of-sample metrics:

gb_final_fit %>% 
  collect_predictions() %>% 
  evaluation_metrics(truth = mwh, estimate = .pred)

Out-of-sample time series:

gb_final_fit %>% 
  extract_workflow() %>% 
  augment(dt_test) %>% 
  transmute(date, prediction = .pred, actual = mwh) %>% 
  pivot_longer(-date) %>% 
  ggplot(aes(date, value, colour = name)) +
  geom_line(alpha = 0.75, size = 0.75) +
  labs(title = "Out-of-sample Time Series",
       y = NULL,
       x = NULL) +
  ggsci::scale_colour_jama() +
  scale_y_continuous(labels = scales::comma_format(suffix = " MWh")) +
  theme_bw() +
  theme(plot.title = element_text(size = 14, face = "bold"),
        plot.subtitle = element_text(face = "italic", size = 12,
                                     colour = "grey50"))

Mean Absolute Error as percentage of average time holdout target value:

mean_abs_error <- gb_final_fit %>%
  extract_workflow() %>% 
  augment(dt_test) %>% 
  mae(mwh, .pred) %>% 
  pull(.estimate)

mean_mwh <- gb_final_fit %>%
  extract_workflow() %>% 
  augment(dt_test) %>% 
  summarise(mean(mwh)) %>% 
  pull()

mean_abs_error/mean_mwh

## [1] 0.01942897

Deviation by date:

gb_final_fit %>% 
  extract_workflow() %>% 
  augment(dt_test) %>% 
  mutate(delta = .pred/mwh-1,
         exceedance = ifelse(abs(delta) > 0.05, 1, 0) %>% as.factor()) %>% 
  ggplot(aes(date, delta, colour = exceedance)) +
  geom_point() +
  ggsci::scale_colour_jama()

gb_final_fit %>% 
  extract_workflow() %>% 
  augment(dt_test) %>% 
  mutate(delta = .pred/mwh-1) %>% 
  transmute(day = yday(date), delta) %>% 
  acf()

Which days still show high deviation? Find out what happened on these days in 2019

gb_final_fit %>% 
  extract_workflow() %>% 
  augment(dt_test) %>% 
  mutate(delta = .pred/mwh-1,
         exceedance = ifelse(abs(delta) > 0.05, 1, 0) %>% as.factor()) %>% 
  filter(exceedance == 1) %>% 
  select(date, exceedance) %>% 
  print(n = 100)

## # A tibble: 15 × 2
##    date       exceedance
##    <date>     <fct>     
##  1 2019-03-29 1         
##  2 2019-04-06 1         
##  3 2019-04-07 1         
##  4 2019-04-11 1         
##  5 2019-04-15 1         
##  6 2019-04-19 1         
##  7 2019-04-22 1         
##  8 2019-05-20 1         
##  9 2019-05-21 1         
## 10 2019-05-31 1         
## 11 2019-06-10 1         
## 12 2019-06-20 1         
## 13 2019-08-01 1         
## 14 2019-12-23 1         
## 15 2019-12-26 1

The tuning managed to improve metrics.

gb_final_fit %>%
  extract_workflow() %>% 
  augment(dt_test) %>% 
  mutate(lower_conf = .pred*(1-mean_abs_error/mean_mwh),
         higher_conf = .pred*(1+mean_abs_error/mean_mwh)) %>%  
  ggplot() +
  geom_line(aes(date, .pred), colour = "midnightblue", alpha = 0.5) +
  geom_line(aes(date, mwh), colour = "firebrick") +
  geom_ribbon(aes(x = date, ymin = lower_conf, ymax = higher_conf), 
              fill = "midnightblue", alpha = 0.25) +
  theme_light()

I hope this post has been interesting to you. In case of constructive feedback or if you want to exchange about this or a related topic, feel free to reach out.

Thank you for reading.

 

---

A work by Mathias Steilen